Since lifting method was first introduced, the lifting method has become a powerful method to perform wavelet transforms and many other orthogonal transforms. Especially for integer-to-integer wavelet transforms, the lifting method is an indispensable tool for lossless image compression. Prior work has shown that the number of lifting steps can have an impact on the transform performance. The fidelity of integer-to-integer transforms depends entirely on how well the integer-to-integer transforms approximate original wavelet transforms. The predominant source of errors in the lifting steps of the lifting method is due to the rounding off of the intermediate real result to an integer at each lifting step. Hence, a wavelet transform with a large number of lifting steps would automatically increase the approximation error. In the case of lossy compression, where a finite-precision arithmetic processor is used, the approximation error is less important because the approximation error is usually masked out by the transform coefficient quantization error. However, in the case of lossless compression, the compression performance is certainly affected by the approximation error. Consequently, the number of lifting steps in a wavelet transform is a major concern.
In conventional prior lifting methods, the arithmetic computation of a single-level and two-channel wavelet transform involves two basic filter convolutions, each of which is followed by a down-sampling operation as shown in filter equations.
            s      ⁡              (        n        )              =                  ∑                  i          =          0                          M          -          1                    ⁢                        h          ⁡                      (            i            )                          ·                  x          ⁡                      (                                          2                ⁢                n                            -              i                        )                                          d      ⁡              (        n        )              =                  ∑                  i          =          0                          M          -          1                    ⁢                        g          ⁡                      (            i            )                          ·                  x          ⁡                      (                                          2                ⁢                n                            -              i                        )                              
In the filter equations, the term x is an input signal, s(n) is the down-sampled output from a lowpass filter h, and d(n) is a down-sampled output from the highpass filter g. For finite impulse response (FIR) filters, lowpass and highpass filters may be split, into even and odd parts defined by h(z)=he(z2)+z−1ho(z2) and g(z)=ge(z2)+z−1go(z2) split equations where subscript e denotes an even part and subscript o denotes an odd part. A P(z) polyphase matrix is composed of even and odd parts of the filter impulse response functions, and is a convenient way to express the wavelet transform operation. Thus, the wavelet transform can now be written as a polyphase wavelet transform matrix equation. In the polyphase wavelet transform matrix equation, P(z) is the polyphase matrix.
            [                                                  s              ⁡                              (                z                )                                                                                        d              ⁡                              (                z                )                                                        ]        =                  P        ⁡                  (          z          )                    ⁡              [                                                                              x                  e                                ⁡                                  (                  z                  )                                                                                                                          z                                      -                    1                                                  ⁢                                                      x                    o                                    ⁡                                      (                    z                    )                                                                                      ]                        P      ⁡              (        z        )              =          [                                                                  h                e                            ⁡                              (                z                )                                                                                        h                o                            ⁡                              (                z                )                                                                                                        g                e                            ⁡                              (                z                )                                                                                        g                o                            ⁡                              (                z                )                                                        ]      
Initially, the conventional lifting method was introduced to improve the performance of wavelet transforms. In a two-channel, single-level wavelet transform, the input signal in each channel is first filtered by convolution and then downsampled by a factor of two. That is, half of the filtered samples are not used and only the other half is used to compute the wavelet transform. The use of only half of the filtered samples is inefficient, and, therefore, downsampling is performed before the filtering in order to reduce the number of operations. The reduction in the number of operations is the basic idea behind the formulation of polyphase matrices. In the polyphase structure of the wavelet transform, the input signal is first split into even and odd parts before entering into the polyphase filters for convolution. By factoring the polyphase matrix into a product of many 2×2 elementary matrices, the lifting method may reduce the convolution operations to an algorithm involving simple multiplication and addition. Many orthogonal and biorthogonal wavelet transforms may be constructed by using some number of iterations of the lifting steps. In particular, Daubechies and Sweldens have presented a systematic analysis of factoring wavelet transforms into lifting steps. The lifting factorization asymptotically reduces the computational complexity of the wavelet transform by a factor of 2. Importantly, the lifting method is the only convenient method that yields integer wavelet coefficients for integer input. The conventional lifting method factorizes the polyphase matrix into a product of 2×2 elementary matrices and a scaling matrix at the end in the form as given by a polyphase factoring equation.
      P    ⁡          (      z      )        =            [                                    c                                0                                                0                                              1              c                                          ]        ⁢                  ∏                  i          =          1                m            ⁢                          ⁢                        [                                                    1                                                                                  t                    i                                    ⁡                                      (                    z                    )                                                                                                      0                                            1                                              ]                ⁡                  [                                                    1                                            0                                                                                                          p                    i                                    ⁡                                      (                    z                    )                                                                              1                                              ]                    
The scalar factor c in the polyphase factoring equation is a normalization factor, usually a floating-point number. The matrix diagonal (c,1/c) in the same equation is called the scaling matrix. Each of the other elementary matrices represents a prediction or an update step. The prediction step elementary matrix
  [                    1                    0                            P                    1              ]is of the form having 1,0 in the first row and p,1 in the second row, and for an update step elementary matrix having 1,t in the first row and 0,1 in the second row, where p and t are polynomials in z with z being the sample index. For most of the cases, the polynomials are zero or first-degree polynomials. The number of elementary matrices for prediction or update is the number of lifting steps or simply as steps. According to the conventional lifting method, the flow of a wavelet transform can thus be stated as splitting the input signal x(i) into even numbered samples x(2i) with s(i)←x(2i) and odd numbered samples x(2i+1) with d(i)←x(2i+1), lifting by the prediction step d(i)←d(i)+Q1(si), lifting by the update step s(i)←s(i)+Q2(di), and normalizing by si←csi and di←di/c, using the scalar c, with Q1(si) and Q2(di) being linear polynomials of si and di, respectively. The flow of an inverse wavelet transform follows the lifting steps in reverse direction. However, even when the input data are integers such as image pixels, the output data from the lifting algorithm are real numbers because of the sum of products in the prediction and update steps, and the final normalization operations. For an integer-to-integer wavelet transform, which is required for lossless data compression, the rounding of the sum of products into an integer in each of the prediction and update steps, but not the scaling factors in the normalization step. Therefore, there is a need to factor the scaling matrix into elementary lifting matrices too. In the prior methods, the scaling matrix is factored into four lifting steps as given by one of the following first, second and third scaling factorization matrices.
      [                            c                          0                                      0                                      1            c                                ]    =                                                        [                                                                    1                                                                              c                      -                      1                                                                                                            0                                                        1                                                              ]                        ⁡                          [                                                                    1                                                        0                                                                                        1                                                        1                                                              ]                                ⁡                      [                                                            1                                                                                            (                                              1                        /                        c                                            )                                        -                    1                                                                                                0                                                  1                                                      ]                          ⁡                  [                                                    1                                            0                                                                                      -                  c                                                            1                                              ]                    ⁢                          [                                    c                                0                                                0                                              1              c                                          ]        =                                                                      [                                                                            1                                                              0                                                                                                                                                    -                          1                                                /                        c                                                                                    1                                                                      ]                            ⁡                              [                                                                            1                                                                                      c                        -                        1                                                                                                                        0                                                              1                                                                      ]                                      ⁡                          [                                                                    1                                                        0                                                                                        1                                                        1                                                              ]                                ⁡                      [                                                            1                                                                                            (                                              1                        -                        c                                            )                                        /                    c                                                                                                0                                                  1                                                      ]                          ⁢                                  [                                            c                                      0                                                          0                                                      1                c                                                    ]            =                                                  [                                                                    1                                                        0                                                                                                              -                      1                                                                            1                                                              ]                        ⁡                          [                                                                    1                                                                              1                      -                                              1                        /                        c                                                                                                                                  0                                                        1                                                              ]                                ⁡                      [                                                            1                                                  0                                                                              c                                                  1                                                      ]                          ⁡                  [                                                    1                                                                                  1                    /                                          c                      2                                                        -                                      1                    /                    c                                                                                                      0                                            1                                              ]                    
Merging the first lifting step of the factorization matrices in scaling matrix with the last lifting step from the rest of wavelet factorization can be accomplished with only three extra steps that are needed to avoid scaling. For example, a simple class of wavelet transforms, such as, (2,2), (2,4), (4,2), (6,2), Haar, and (4,4) biorthogonal wavelet transforms may be factored into polyphase factoring equations. The polyphase matrix in the conventional lifting method of such wavelet transforms consists of three elementary matrix factors including one prediction step, one update step, and one scaling matrix.
      P    ⁡          (      z      )        =                    [                                            c                                      0                                                          0                                                      1                c                                                    ]            ⁡              [                                            1                                                      Q                2                                                                        0                                      1                                      ]              ⁡          [                                    1                                0                                                              Q              1                                            1                              ]      
In the polyphase matrix of (2,2) biorthogonal wavelet transform, for example, c=√2, Q1=−(1+z−1)/2, and Q2=(1+z)/4. The sample index z−1 denotes the next sample and z, the immediately preceding sample. Therefore, for (2,2) biorthogonal wavelet transform, Q1(si)=−{(s(i)+s(i+1)}/2 and Q2(di)={d(i)+d(i−1)}/4. Neglecting the scaling matrix, the integer-to-integer lifting process of the (2,2) biorthogonal wavelet transform comprises of one prediction and one update lifting steps. In the initial step the lifting process is to split the input samples x(i) into even numbered and odd numbered samples as denoted by s(0)(i) and d(0)(i) respectively as s(0)(i)=x(2i) and d(0)(i)=x(2i+1). The following prediction step is to give a new odd numbered sample as d(1)(i)=d(0)(i)−[(s(0)(i)+s(0)(i+1))/2]. The next update step is to calculate an updated even numbered sample as s(1)(i)=s(0)(i)+[(d(1)(i)+d(1)(i−1))/4]. The square brackets [ . . . ] denotes rounding-off to the nearest integers. The value within the [ ] at the prediction step is a linear interpolation of two adjacent s(0)(i) and s(0)(i+1) values. Thus the d(1)(i) is the error between the odd numbered sample value, d(0)(i), and the linearly interpolated value of two adjacent even numbered samples.
In terms of the original input samples x(i), the d(1)(i) value is given by d(1)(i)=[(−x(2i)+2x(2i+1)−x(2i+2))/2]. Also, in terms of the original samples, the updated value at the update step is obtained as s(1)(i)=[(−x(2i−2)+2x(2i−1)+6x(2i)+2x(2i+1)−x(2i+2))/8]. These two samples, d(1)(i) and s(1)(i) are just the highpass and the lowpass outputs of the 5/3 biorthogonal wavelet transform used by the international image compression standard, JPEG 2000 for lossless image data compression.
For integer-to-integer transform, which is required in lossless data compression, the scaling matrix in the polyphase matrix can be replaced by the one of the factorization matrices, for example, the third factorization matrix to derive an expanded P(z) matrix.
      P    ⁡          (      z      )        =                                          [                                                            1                                                  0                                                                                                  -                    1                                                                    1                                                      ]                    ⁡                      [                                                            1                                                  ck                                                                              0                                                  1                                                      ]                          ⁡                  [                                                    1                                            0                                                                    c                                            1                                              ]                    ⁡              [                                            1                                                                        Q                  2                                -                k                                                                        0                                      1                                      ]              ⁡          [                                    1                                0                                                              Q              1                                            1                              ]      
In the expanded P(z) matrix, k=(c−1)/c2. Therefore, the conventional lifting method disadvantageously requires five lifting steps for (2,2), (2,4), (4,2), (6,2), Haar, and (4,4) biorthogonal wavelet transforms.
The polyphase matrix of another class of wavelet transforms, such as, the fourth order Daubechies orthogonal wavelet (D4), the Sequential and Predictive (S+P), and the (1,1+1), TS, (2,10), and (2+2,2) biorthogonal wavelet transforms. The wavelet transforms consists of four elementary matrix factors including one first prediction step, one update step, one second prediction step, and one scaling matrix.
      P    ⁡          (      z      )        =                              [                                                    c                                            0                                                                    0                                                              1                  c                                                              ]                ⁡                  [                                                    1                                            0                                                                                      Q                  3                                                            1                                              ]                    ⁡              [                                            1                                                      Q                2                                                                        0                                      1                                      ]              ⁡          [                                    1                                0                                                              Q              1                                            1                              ]      
For integer-to-integer transform, which is required in lossless data compression, the scaling matrix in the polyphase matrix is replaced by one of the scaling factorization matrices. Therefore, the conventional lifting method disadvantageously requires six lifting steps for the lossless fourth-order Daubechies orthogonal wavelet (D4), S+P, and the (1,1+1), TS, (2,10), and (2+2,2) biorthogonal wavelet transforms.
The polyphase matrix of yet another class of wavelet transforms, such as, the sixth-order Daubechies (D6) orthogonal wavelet transform, and the 9/7 biorthogonal wavelet transform, that consists of five elementary matrix factors. These wavelet transforms include one first prediction step, a first update step, a second prediction step, a second update step, and a scaling matrix. The conventional lifting method disadvantageously requires seven lifting steps for the sixth order Daubechies orthogonal (D6) wavelet transform and the 9/7 biorthogonal wavelet transform.
In a second aspect of the prior art, the performance of using lifting method for lossless data compression has been relatively unpredictable. The reasons are due, firstly, to the overall rounding errors accumulated in many lifting steps that diminish the approximation power of wavelet transforms, and secondly, to the lifting scheme that is not locally adaptive to image properties. The advantage of using a wavelet transform for image data compression is the power of adapting to local properties of the pixels. In remote sensing hyperspectral data, many but not all spectral planes are well correlated. In each spectral plane, the spatial data is composed of patches of relatively smooth areas segmented by edges. The smooth areas can be well compressed by a relatively long wavelet transform with a larger number of vanishing moments. However, for the regions around edges, shorter wavelet transforms are preferable. Despite the fact that the local statistics of both the spectral and spatial data change from pixel to pixel, almost all known image data compression algorithms use only one wavelet transform for the entire dataset. For example, the current international image data compression standard, JPEG 2000, has adopted the (2,2) wavelet transform as the default standard for lossless image data compression for all still images. There is not a single wavelet filter that performs uniformly better than the others. Thus, it would be beneficial to use many types of wavelet filters based on local activities of the image. The selected wavelet transform can thus be best adapted to the content of the image locally.
The lifting method is a fast and powerful tool to implement all wavelet transforms. Especially for integer-to-integer wavelet transforms, the lifting method is an indispensable tool for lossless image compression. Many reversible integer-to-integer wavelet transforms have been evaluated for lossless image compression, but no adaptive lifting scheme has been mentioned. Edge adaptive lifting methods have switched the orders of prediction and update steps when applied to lossy image data compression. An adaptive lifting method has been used for lossless image data compression. The prediction step is the median of six predictors that were parts of the prediction filters of the biorthogonal wavelet transforms (2,2), (4,2), and (6,2), but poor compression ratio performance for the same image quality is disadvantageously expected from poor prediction image edges.
The current international image data compression standard, JPEG 2000, has adopted the (2,2) biorthogonal wavelet transform as the default standard of lossless image data compression for all still images. The primary reason for such a decision is that the (2,2) biorthogonal wavelet filter is relatively short and can closely follow the variation of local image contents. The lifting scheme of the (2,2) wavelet transform consists of one prediction step and one update step. At the prediction step, a linear interpolation of two adjacent even-numbered samples is used to get an estimate of the odd-numbered sample in the middle. FIG. 1 depicts conventional linear interpolations of four possible scenarios. In FIGS. 1A and 1B the slopes at the two adjacent even-numbered samples have the same sign, whereas in FIGS. 1C and 1D the two adjacent even-numbered samples have the opposite sign. As shown, the linear interpolation is more accurate in the former where the two slopes have the same sign than the latter where the two slopes have the opposite sign. In the latter case, a cubic interpolation would be more accurate. The estimated odd-numbered sample due to linear interpolation of two adjacent even-numbered samples is given by linear=(si+si+1)/2. The estimated odd-numbered sample due to a four-point cubic interpolation of four adjacent even-numbered samples is given by cubic=(−si−1+9si+9si+1−si+2)/16. The first order difference between the cubic and linear interpolated values is Δ(i)=(si−1−si−si+1+si+2)/16.
Neglecting the scaling factor, the conventional (2+2,2) integer-to-integer lifting wavelet transform comprises of one first prediction, one first update, and a second prediction lifting steps. In the initial step, the lifting process is to split the input samples x(i) into even numbered and odd numbered samples as denoted by s(0)(i) and d(0)(i) respectively as s(0)(i)=x(2i) and d(0)(i)=x(2i+1). The following first prediction step is to give a new odd numbered sample as d(1)(i)=d(0)(i)−[(s(0)(i)+s(0)(i+1))/2]. The next update step is to calculate an updated even numbered sample as s(1)(i)=s(0)(i)+[(d(1)(i)+d(1)(i−1))/4]. The following second prediction step is to give a new odd numbered sample as d(2)(i)=d(1)(i)−[(s(1)(i)−s(1)(i−1)+s(1)(i+1)−s(1)(i+2))/16]. The first three lifting steps in (2+2,2) biorthogonal wavelet, split, prediction, and update, are the same as in the conventional integer-to-integer lifting of (2,2) biorthogonal wavelet transform. Thus the predicted value at the first prediction step in the conventional (2,2) and (2+2,2) wavelets is a linear interpolated value as given by linear=(si+si+1)/2. The predicted value at the second prediction step in the conventional (2+2,2) wavelet is the first order difference between the cubic and linear interpolated values, Δ(i). The first order difference is a mathematical term for commonly named difference. In terms of sample index, the operation of Δ(z) is written as Δ=(z−1−z−1+z−2)/16. At the second prediction step in wavelet (2+2,2) the estimated odd-numbered sample is the difference between cubic and linear interpolations. However, the cubic and linear interpolation difference value at the second prediction step is not good enough, because the first order difference does not take into account the sudden changes in the first update result.
In a third aspect of the prior art, the current international standard for still image compression, JPEG 2000, has two completely different default wavelet transforms, one is the 5/3 biorthogonal wavelet transform used for lossless data compression, and the other, CDF 9/7 biorthogonal wavelet transform for lossy data compression. The highpass filter coefficients of the 5/3 biorthogonal wavelet transforms, in terms of rational numbers are given by h1=(−1, 2, −1)/2, and the lowpass filter, h0=(−1, 2, 6, 2, −1)/8. The highpass filter coefficients of the CDF 9/7 biorthogonal wavelet transforms, in terms of rational numbers are given by h1=(1, 0, −9, 16, −9, 0, 1)/16, and the lowpass filter, h0=(1, 0, −8, 16, 46, 16, −8, 0, 1)/64.
Due to the completely different highpass and lowpass filter coefficients in the 5/3 and CDF 9/7 biorthogonal wavelet transforms, the polyphase matrices for the 5/3 and CDF 9/7 biorthogonal wavelet transforms must be factorized into two completely different elementary matrices. Thus, the 5/3 and CDF 9/7 biorthogonal wavelet transforms cannot be implemented using a common lifting step, and hence cannot share processor resources, as such, portion of the lossless lifting method cannot be used for the lossy lifting method. In some downlink applications communicating lossless compression data, the received compressed data must be first disadvantageously decompressed back into the original data set and then lossy compressed in order to generate a lossy compression data set for bandwidth efficient communication of images of a local communications network for broadband distribution. These and other disadvantages are solved or reduced using the invention.